TIME REVERSAL FOR WAVES IN RANDOM MEDIA 
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Abstract. In time reversal acoustics experiments, a signal is emitted from a localized source, 
recorded at an array of receivers-transducers, time reversed, and finally re-emitted into the medium. 
A celebrated feature of time reversal experiments is that the refocusing of the re-emitted signals 
at the location of the initial source is improved when the medium is heterogeneous. Contrary to 
intuition, multiple scattering enhances the spatial resolution of the refocused signal and allows one 
to beat the diffraction limit obtained in homogeneous media. This paper presents a quantitative 
explanation of time reversal and other more general refocusing phenomena for general classical waves 
in heterogeneous media. The theory is based on the asymptotic analysis of the Wigner transform of 
wave fields in the high frequency limit. Numerical experiments complement the theory. 
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1. Introduction. In time reversal experiments, acoustic waves are emitted from 
a localized source, recorded in time by an array of receivers-transducers, time reversed, 
and re-transmitted into the medium, so that the signals recorded first are re-emitted 
last and vice versa 0, [|, |l3|, [l6| [H| |f9). The re-transmitted signal refocuses at the 
location of the original source with a modified shape that depends on the array of 
receivers. The salient feature of these time reversal experiments is that refocusing 
is much better when wave propagation occurs in complicated environments than in 
homogeneous media. Time reversal techniques with improved refocusing in heteroge- 
neous medium have found important applications in medicine, non-destructive testing, 
underwater acoustics, and wireless communications (see the above references). It has 
been also applied to imaging in weakly random media [l^] 
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A schematic description of the time reversal procedure is depicted in Fig. 
Early experiments in time reversal acoustics are described in Q; see also the more 
recent papers [[ll], [l2| |T^]. A very qualitative explanation for the better refocusing 
observed in heterogeneous media is based on multipathing. Since waves can scatter 
off a larger number of heterogeneities, more paths coming from the source reach the 
recording array, thus more is known about the source by the transducers than in a 
homogeneous medium. The heterogeneous medium plays the role of a lens that widens 
the aperture through which the array of receivers sees the source. Refocusing is also 
qualitatively justified by ray theory (geometrical optics). The phase shift caused by 
multiple scattering is exactly compensated when the time reversed signal follows the 
same path back to the source location. This phase cancellation happens only at 
the source location. The phase shift along paths leading to other points in space is 
essentially random. The interference of multiple paths will thus be constructive at 
the source location and destructive anywhere else. This explains why refocusing at 
the source location is improved when the number of scatterers is large. 

As convincing as they are, the above explanations remain qualitative and do not 
allow us to quantify how the refocused signal is modified by the time reversal proce- 
dure. Quantitative justifications require to analyze wave propagation more carefully. 
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Fig. 1.1. The Time Reversal Procedure. Top: Propagation of signal and measurements in time. 
Bottom: Time reversal of recorded signals and back-propagation into the medium. 



The first quantitative description of time reversal was obtained in |j in the framework 
of randomly layered media (see also the recent work 0). That paper provides the 
first mathematical explanation of two of the most prominent features of time reversal: 
heterogeneities improve refocusing and rcfocusing occurs for almost every realization 
of the random medium. The first multi-dimensional quantitative description of time 
reversal was obtained in Q] for the parabolic approximation, i.e., for waves that prop- 
agate in a privileged direction with no backscattering (see also J23| for further analysis 
of time reversal in this regime). That paper shows that the random medium indeed 
plays the role of a lens. The back- propagated signal behaves as if the initial array 
were replaced by another one with a much bigger effective aperture. In a slightly 
different context, a recent paper Q analyzes time reversal in ergodic cavities. There, 
wave mixing is created by reflection at the boundary of a chaotic cavity, which plays 
a similar role to the heterogeneities in a heterogeneous medium. 

This paper generalizes the results of ||] to the case of general classical waves 
propagating in weakly fluctuating random media. The main results are briefly sum- 
marized as follows. We first show that rcfocusing in time reversal experiments may 
be understood in the following three-step more general framework: 

(i) A signal propagating from a localized source is recorded at a single time T > 
by an array of receivers. 

(ii) The recorded signal is processed at the array location. 
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(iii) The processed signal is emitted from the array and propagates in the same 
medium during the same time T. 
The first main result is that the resulting signal will refocus at the location of the 
original source for a large class of waves and a large class of processings. The ex- 
periments described above correspond to the specific processing of acoustic waves in 
which pressure is kept unchanged and the direction of the acoustic field is reversed. 

The second main result is a quantitative description of the re-transmitted signal. 
We show that the re-propagated signal u s (£) at a point £ near the source location 
can be written in the high frequency limit as the following convolution of the original 
source S 

u B (£) = (F*SM). (1.1) 

The kernel F depends on the location of the recording array and on the signal pro- 
cessing. The quality of the refocusing depends on the spatial decay of F. It turns 
out that it can be expressed in terms of the Wigner transform |24|] of two wave fields. 
The decay properties of F depend on the smoothness of the Wigner transform in the 
phase space. The Wigner transform in random media has been extensively studied 
|), 24, 2(|, especially in the high frequency regime, when the wavelength of the initial 



signal is small compared to the distance of propagation. It satisfies a radiative trans- 
port equation, which is used to describe the evolution of the energy density of waves 
in random media [|l7| [2l|. |2^, |2G| . The transport equations possess a smoothing effect 
so that the Wigner distribution becomes less singular in random media, which implies 
a stronger decay of the convolution kernel F and a better refocusing. The diffusion 
approximation to the radiative transport equations provides simple reconstruction 
formulas that can be used to quantify the refocusing quality of the back-propagated 
signal. This construction applies to a large class of classical waves: acoustic, electro- 
magnetic, elastic, and others, and allows for a large class of signal processings at the 
recording array. 

Some results of this paper have been announced in ]IJ . The concept of one-step 
time reversal emerged during early discussions with Knut Solna. We also stress that 
the important property of self-averaging of the time reversed signal (the refocused 
signal is almost independent of the realization of the random medium) is not analyzed 
in this paper. A formal explanation is given in [Q, [2^] in the parabolic approximation. 
Self-averaging for classical waves will be addressed elsewhere. 

This paper is organized as follows. Section ^| recalls the classical setting of time 
reversal and introduces one-step time reversal. The re-transmitted signal and its rela- 
tion to the Wigner transform are analyzed in section ^. A quantitative description of 
acoustic wave refocusing in weakly fluctuating random media is obtained by asymp- 
totic analysis; see equations (3.42) and ( |3.43| ) for an explicit expression in the diffusion 



approximation. Section |] generalizes the results in two ways. First, a more general 
signal processing at the recording array is allowed, such as recording only the pressure 
field of acoustic waves and not the velocity field. Second, the re-transmission scheme 
is applied to more general waves and the role of polarization and mode coupling is 
explained. 

We would like to thank Knut Solna for fruitful discussions during the preparation 
of this work. We are indebted to George Papanicolaou for his contributions to the 
analysis of time reversal, which lie at the core of this paper. This work would also not 
have been possible without the numerous exchanges we benefited from at the Stanford 
MGSS summer school. 
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2. Classical Time Reversal and One-Step Time Reversal. Propagation 
of acoustic waves is described by a system of equations for the pressure p(t, x) and 
acoustic velocity v(i, x): 

<9v 

p(x)— +V P = Q (2.1) 
K (x)^+V-v = 0, 

with suitable initial conditions and where p(x) and /c(x) are density and compress- 
ibility of the underlying medium, respectively. These equations can be recast as the 
following linear hyperbolic system 

, , . du ■ du , 

A ^ai +DJ W = ' xeM (2 ' 2) 

with the vector u = (v,p) € C 4 . The matrix A = Diag(p, p, p, k) is positive definite. 
The 4 x 4 matrices Z3 J , j = 1, 2, 3, are symmetric and given by £>4„ = <5 m 4<5nj+<5n4<i>mj- 
We use the Einstein convention of summation over repeated indices. 

The time reversal experiments in |Q consist of two steps. First, the direct problem 

u(0,x) = S(x) 

with a localized source S centered at a point xo is solved. The signal is recorded 
during the period of time < t < T by an array of receivers located at f2 C R 3 . 
Second, the signal is time reversed and re-emitted into the medium. Time reversal 
is described by multiplying u = (v,p) by the matrix T = Diag(— 1, —1, —1, 1). The 
back-propagated signal solves 

^+A- 1 (x) J D^ = lR(2T-t,x), T<i<2T (2.4) 
u(T,x) = 

with the source term 

R(t,x)=ru(t,x)x(x). (2.5) 

The function x(x) is either the characteristic function of the set where the recording 
array is located, or some other function that allows for possibly space-dependent 
amplification of the re-transmitted signal. 

The back-propagated signal is then given by u(2T, x). We can decompose it as 



u(2T» = i f ds w(s,x;s) 
1 Jo 



(2.6) 



where the vector- valued function w(f, x; s) solves the initial value problem 
w(0,x; s) = R(s,x). 
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Wc deduce from (2.6) that it is sufficient to analyze the refocusing properties of 
w(s,x;s) for < s < T to obtain those of u(2T, x). For a fixed value of s, we call 
the construction of w(s, x; s) one-step time reversal. 

We define one-step time reversal more generally as follows. The direct problem 
( |2.3| ) is solved until time t = T to yield u(T _ ,x). At time T, the signal is recorded 
and processed. The processing is modeled by an amplification function x( x )j a blur- 
ring kernel /(x), and a (possibly spatially varying) time reversal matrix T. After 
processing, we have 



u(T H 



r(/*(xu))(T-,x)x(x). 



The processed signal then propagates during the same time T: 



du du 



u(T+,x) 



0, T < t < 2T 
r(/*( X u))(T- x)x(x). 



(2.7) 



(2.8) 



The main question is whether u(2T, x) refocuses at the location of the original source 
S(x) and how the original signal has been modified by the time reversal procedure. 
Notice that in the case of full (f2 = R 3 ) and exact (/(x) = <5(x)) measurements 
with r = Diag(— 1, —1, —1, 1), the time-reversibility of first-order hyperbolic systems 
implies that u(2T, x) = TS(x), which corresponds to exact refocusing. When only 
partial measurements are available we shall see in the following sections that u(2T, x) 
is closer to TS(x) when propagation occurs in a heterogeneous medium than in a 
homogeneous medium. 

The pressure field p(t, x) satisfies the following scalar wave equation 



d 2 p 



1 



-V • 



k(x) Vp( x ) 



1 



Vp = 0. 



(2.9) 



A schematic description of the one-step procedure for the wave equation is presented 
in Fig. |2.l[ This is the equation solved in the numerical experiments presented 
in this paper. The details of the numerical setting are described in the appendix. A 
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Fig. 2.1. The One-Step Time Reversal Procedure. Here, pt denotes 



dp 

at' 



numerical experiment for the one-step time reversal procedure is shown in Fig. 2^2. In 
the numerical simulations, there is no blurring, /(x) = 5(x), and the array of receivers 
is the domain SI = (—1/6, 1/6) 2 (x(x) is the characteristic function of f2). Note that 
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Fig. 2.2. Numerical experiment of the one-step time reversal procedure. Top Left: initial con- 
dition p(0,x), a peaked Gaussian of maximal amplitude equal to 1. Top Right: forward solution 
p(T~,x), of maximal amplitude 0.04. Bottom Right: recorded solution p(T+ , x) , of maximal am- 
plitude 0.015 on the domain fl = (— 1/6, 1/6) 2 . Bottom Left: back-propagated solution p(2T, x), o/ 
maximal amplitude 0.07. 



the truncated signal does not retain any information about the ballistic part (the part 
that propagates without scattering with the underlying medium). In homogeneous 
medium, the truncated signal would then be identically zero and no refocusing would 
be observed. The interesting aspect of time reversal is that a coherent signal emerges 
at time 2T out of a signal at time T + that seems to have no useful information. 

3. Theory of Time Reversal in Random Media. Our objective is now to 
present a theory that explains in a quantitative manner the refocusing properties 
described in the preceding sections. We consider here the one-step time reversal for 
acoustic wave. Generalizations to other types of waves and more general processings 



in ( |2.8| ) are given in section | 



3.1. Refocused Signal. We recall that the one-step time reversal procedure 
consists of letting an initial pulse S(x) propagate according to ([2.3|) until time T, 



u(T-,x)=/ G(T,x;z)S(z)dz, 

Jk 3 
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where G(T, x; z) is the Green's matrix solution of 

gG(tx;y) .9Gfex;y) = 

G(0,x;y)=/S(x-y). 

At time T, the "intelligent" array reverses the signal. For acoustic pulses, this means 
keeping pressure unchanged and reversing the sign of the velocity field. The array 
of receivers is located in O C R 3 . The amplification function x( x ) is an arbitrary 
bounded function supported in £1, such as its characteristic function (x( x ) = 1 for 
x G il and x( x ) = otherwise) when all transducers have the same amplification 
factor. We also allow for some blurring of the recorded data modeled by a convolution 
with a function /(x). The case /(x) = <5(x) corresponds to exact measurements. 
Finally, the signal is time reversed, that is, the direction of the acoustic velocity is 



reversed. Here, the operator T in (2.7) is simply multiplication by the matrix 

r = Diag(-l, -1,-1,1). (3.2) 
The signal at time T + after time reversal takes then the form 

u(T+,x)=/ rG(T,y';z) X (x) X (y')/(x-y')S(z)dzdy'. (3.3) 



The last step (2.S) consists of letting the time reversed field propagate through 
the random medium until time 2T. To compare this signal with the initial pulse S, 
we need to reverse the acoustic velocity once again, and define 

u B (x) = ru(2T,x) = f rG(r,x; y )rG(r, y ';z)x(y)x(y')/(y-y')s(z)dyd y 'dz. 

(3.4) 

The time reversibility of first-order hyperbolic systems implies that u s (x) = 
S(x) when £1 = R 3 , \ = 1, and /(x) = <5(x), that is, when full and non-distorted 
measurements are available. It remains to understand which features of S are retained 
by u B (x) when only partial measurement is available. 

3.2. Localized Source and Scaling. We consider an asymptotic solution of the 



time reversal problem (2.3), ( |2.8| ) when the support A of the initial pulse S(x) is much 
smaller than the distance L of propagation between the source and the recording array: 
£ = X/L <C 1. We also take the size a of the array comparable to L: a/L = 0(1). 
We assume that the time T between the emission of the original signal and recording 
is of order L/co, where cq is a typical speed of propagation of the acoustic wave. We 
consequently consider the initial pulse to be of the form 

u(0,x) = S(^^) 

e 

in non-dimensionalized variables x' = x/L and t' = t/(L/co). We drop primes to 
simplify notation. Here Xo is the location of the source. The transducers obviously 
have to be capable of capturing signals of frequency and blurring should happen 
on the scale of the source, so we replace /(x) by e _3 /(e _1 x). Finally, we are interested 



8 



G. BAL and L. RYZHIK 



in the refocusing properties of u (x) in the vicinity of x . We therefore introduce the 



scaling x = x + e£. With these changes of variables, expression (3.4) is recast as 

u B (£;x ) =ru(2T,xo + e£) (3.5) 
= / rG(r,x +^;y)rG(T,y';x +ez) X (y,y')S( Z )ciydy'dz, 



where 



x(y,y') = x(y)x(y')fO L -^)- (3.6) 



In the sequel we will also allow the medium to vary on a scale comparable to the 
source scale e. Thus the Green's function G and the matrix A depend on e. We 
do not make this dependence explicit to simplify notation. We are interested in the 
limit of u B (£; xo) as e — ► 0. The scaling considered here is well adapted both to the 
physical experiments in and the numerical experiments in Fig. |2.2| . 

3.3. Adjoint Green's Function. The analysis of the re-propagated signal relies 
on the study of the two point correlation at nearby points of the Green's matrix in 



(3.5). There are two undesirable features in (3.5). First, the two nearby points Xn + e£ 
and xo+ez are terminal and initial points in their respective Green's matrices. Second, 
one would like the matrix T between the two Green's matrices to be outside of their 
product. However, T and G do not commute. For these reasons, we introduce the 
adjoint Green's matrix, solution of 

dG^-y) aG.(t,x;y) 

dt v ; dxi (3-7) 

G <! (0,x;y)=^- 1 (x)<5(x-y). 

We now prove that 

G4t,x;y)=rG(i,y;x)A- 1 (x)r. (3.8) 
Note that for all initial data S(x), the solution u(t, x) of (2.3) satisfies 

u(f,x)=/ G(t - s, x; y)u(s, y)dy 

for all < s < t < T since the coefficients in (2.3) are time-independent. Differenti- 



ating the above with respect to s and using (2.3) yields 
Upon integrating by parts and letting s = 0, we get 

Since the above relation holds for all test functions S(y), we deduce that 

fez)-|_ [G(f , x;y) ^ 1(y)DJ] .„. 



(3.9) 
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Interchanging x and y in the above equation and multiplying it on the left and the 
right by T, we obtain that 

^ [TG(t, y; x)^ 1 (x)] A(x)T - A [TG(t, y; x)^ 1 (x)] EfiT = 0. (3.10) 
We remark that 

TD j = -D j T and TA(x) = A(x)T, (3.11) 

so that 

^ [rGfryrfA-^r] A(x) + A [rG(f, y; x)^ 1 (x)r] D J = 

with rG(0,y;x) J 4- 1 (x)r = A- 1 (x)S{x - y). Thus (|J) follows from the uniqueness 
of the solution to the above hyperbolic system with given initial conditions. We can 
now recast ( |3.5| ) as 

u B (£;x ) = / rG(T,xo+££;y)G*(T,xo + £z;y')r 

Jk 9 _ , (3.12) 

x X (y) X (y , )/(^ L )^(xo + ez)S(z)dydy'dz. 

One further simplifies ( |3.12| ) with the help of the auxiliary matrix- valued functions 
Q(t, x;q) and Q*(i,x, q) defined by 

Q(T,x;q) = f G(T,x;y)x(y)e^/ £ dy, 

Jr 3 (3.13) 

Q*(7>;q) = / G,(T,x;y) X (y)e- iqy / £ dy. 
it 3 

They solve the hyperbolic equations ( |2.3| ) and ( |3.7| ) with initial conditions given by 
Q(0,x;q) = x(x)e lc > x / 6 / and Q*(0,x;q) = A~i(x)x(x)e~ ?q ' x/e , respectively. Thus 
( p. 12 ) becomes 



u s (£; x )= / rQ(T, x + e€; q)Q.(T, x + ez; q)TA(x Q + ez)S(z)/(q) 



dqdz 

(2^)3' 

(3.14) 



where /(q) = L 3 e tqx /(x)<ix is the Fourier transform of /(x). 

3.4. Wigner Transform. The back-propagated signal in ( 3.14 ) now has the 
suitable form to be analyzed in the Wigner transform formalism fl4j, |24f| . We define 

W e (t,x,k)= / /(q)[/ £ (i,x,k;q)dq, (3.15) 

Jr 3 



where 



U e (t,x,k;q)= [ e ' :k -"Q(f,x^;q)Q,(f,x+^;q)A. ( 3 .i 6 ) 
J R3 2 2 (2tt) j 

Taking inverse Fourier transform we verify that 

0(*,x;q)Q,(t,y;q)= / e -* (y- x )/-[/ e (t, ^±Z, k; q)dk, 

Jm 3 * 
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hence 

dzdk 

W) 



u B (|;x )= f e ik -(«-)r^(T,x + £^,k)rA(x + e z)S(z)^ 



We have thus reduced the analysis of u(£; x ) as e — » to that of the asymptotic 
properties of the Wigner transform W e . The Wigner transform has been used exten- 
sively in the study of wave propagation in random media, especially in the derivation 
of radiative transport equations modeling the propagation of high frequency waves. 
We refer to [O, |2lL e3. Note that in the usual definition of the Wigner transform, 



one has the adjoint matrix Q* in place of Q* in ( 3.16 ). This difference is not essential 
since and Q* satisfy the same evolution equation, though with different initial 
data. 



The main reason for using the Wigner transform in (3.17) is that W e has a weak 
limit W as e — > 0. Its existence follows from simple a priori bounds for W E (t, x, k). 
Let us introduce the space A of matrix-valued functions 0(x, k) bounded in the norm 
|| • \\ A defined by 

||0IU= / sup||0(x,y)||dy, where 0(x,y) = / e~ 4k ' y 0(x, k)dk. 

We denote by A' its dual space, which is a space of distributions large enough to 
contain matrix-valued bounded measures, for instance. We then have the following 
result: 

Lemma 3.1. Let x(x) € L 2 (R 3 ) and /(q) £ X 1 (R 3 ). Then there is a constant 
C > independent of e > and t £ [0, oo) such that for all t S [0, oo), we have 

iiw e (t > x,k)m,<c. 

The proof of this lemma is essentially contained in |l4|, pjj ; sec also Q|. One may 
actually get L 2 -bounds for W e in our setting because of the regularizing effect of / in 



(3.15) but this is not essential for the purposes of this paper. We therefore obtain the 
existence of a subsequence £k — > such that W Ek converges weakly to a distribution 
W G ^4'. Moreover, an easy calculation shows that at time t = 0, we have 

t¥(0,x o ,k) - | X (x )| 2 A - 1 (x )/>). (3.18) 

Here, Aq = A when A is independent of e, and Ao = lim A £ if we assume that the 

family of matrices A £ ('x.) is uniformly bounded and continuous with the limit Aq in 
C(M. d ). These assumptions on A e are sufficient to deal with the radiative transport 



regime we will consider in section 3.7. Under the same assumptions on A e , we have 
the following result. 



Proposition 3.2. The back-propagated signal u ij (^;x ) given by (3.11) con- 
verges weakly in S' (M 3 xl 3 ) as e — > to the limit 

u s (£;xo)=/ e tk («- z 'r^(T,x ,k)rA (x )S(z)^. (3.19) 

The proof of this proposition is based on taking the duality product of u s (£;xo) 
with a vector-valued test function </>(£; xo) in 6>(R 3 x R 3 ). After a change of variables 
we obtain (u B ,<p) = (W e ,Z £ ). Here the duality product for matrices is given by the 
trace (A, B) = J2i,k( A ik, B ik), and 

Z B (xo,k)=/ e 4k -( z -«)r^,x -e^)S*(z)A e (x 0+£ ^)r-^|. (3.20) 

7 R 6 2 2 (27T)' 3 
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Defining Z as th e limit of Z e as e — > by replacing formally e by in the above 
expression, ( 3.19 ) follows from showing that \\Z E — Z\\a — ► as e — ► 0. This is 
straightforward and we omit the details. 

The above proposition tells us how to reconstruct the back-propagated solution 
in the high frequency limit from the limit Wigner matrix W. Notice that we have 
made almost no assumptions on the medium described by the matrix A e (x). At this 
level, the medium can be either homogeneous or heterogeneous. Without any further 
assumptions, we can also obtain some information about the matrix W. Let us define 
the dispersion matrix for the system (|1|) as [E4| 



L(x,k) = A 1 (x)k j D : >. 



(3.21) 



It is given explicitly by 



L(x, k) 



fcx/p(x)\ 
fc 2 /p(x) 
fe/p(x) 
/ 









\fci/«(x) k 2 /n(x) fc 3 //c(x) 

The matrix L has a double eigenvalue loq = and two simple eigenvalues w± (x, k) = 
±c(x)|k|, where c(x) = 1/y p(x)k(x) is the speed of sound. The eigenvalues o->± are 
associated with eigenvectors b± (x, k) and the eigenvalue too = is associated with 
the eigenvectors b,(x, k), j = 1,2. They are given by 



/ 



b±(x,k) 



±- 



V>(x) 



\ 



bj(x, k) 



z J (k) 
o 



(3.22) 



where k = k/|k| and z 1 (k) and z 2 (k) are chosen so that the triple (k, z 1 (k), z 2 (k)) 
forms an orthonormal basis. The eigenvectors are normalized so that 



(A (x)b 3 (x,k) • b fc (x,k)) = S jk , 



(3.23) 



for all j, k G J = {+, — , 1, 2}. The space of 4 x 4 matrices is clearly spanned by the 
basis hj (g) bfc. We then have the following result: 

Proposition 3.3. There exist scalar distributions a± and a nn , m,n = 1,2 so 
that the limit Wigner distribution matrix can be decomposed as 

2 



W(t,x,k)= J2 ^ m (t,x,k)b J -(x,k)(g)b m (x,k) 

j,m=l 

+ o+(t,x, k)b+(x, k) ®b+(x, k) + a_ (t, x, k)b_ (x, k) ®b_(x,k). 



(3.24) 



The main result of this proposition is that the cross terms b^ ® hk with uij ^ Wk do 
not contribute to the limit W. The proof of this proposition can be found in |l4j and 
a formal derivation in p4) . 

The initial conditions for the amplitudes aj are calculated using the identity 



Aq 1 (x) = h 3 ( x > k ) ® b J ( x > k ) • 



Then (|3.18[) implies that aj 2 (0,x,k) = 0^(0, x,k) = and 

a^'(0,x,k) = a±(0,x,k) = |x(x)| 2 /(k), j = 1,2. 



(3.25) 
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3.5 . M ode Decomposition and Refocusing. We can use the above result to 
recast ( [3.19 ) as 



u B (£;x ) = OF(T,.;x )*S)(0, 



(3.26) 



where 



dk 

(2^)3 



V f e i;k -«a™ n (^x ;k)rb m (x ,k)®b„(x ,k)A (x )r 

f dk 

/ e tk -« a+ (T,x ;k)rb + (xo,k)®b + (x ,k)A (x )r— - (3.27) 

f dk 

/ e« k « a _(T, x ; k)rb_ (x , k) ® b_(x , k)A (x )r-— 3 . 



This expression can be used to assess the quality of the refocusing. When F(T, £; x ) 
has a narrow support in £, refocusing is good. When its support in £ grows larger, its 
quality degrades. The spatial decay of the kernel F(t, £; xo) in £ is directly related to 
the smoothness in k of its Fourier transform in £: 



dk 

(2^ 



F(T,k;x ) = J2 ao nn (^x ;k)rb m (x ,k)®b„(x ,k)A (x )r 

m,n— 1 

+r [ a+ (T, x ; k)b+ (x , k) ® b+ (x , k) + a_ (T, x ; k)b (x , k) ® b_ (x , k)] A (x )r. 

Namely, for F to decay in £, one needs F(k) to be smooth in k. However, the 
eigenvectors hj are singular at k = as can be seen from the explicit expressions 
(3.22). Therefore, F a priori is not smooth at k = 0. This means that in order 
to obtain good refocusing one needs the original signal to have no low frequencies: 
S(k) = near k = 0. Low frequencies in the initial data will not refocus well. 

We can further simplify (^26|)-(^27|) is we assume that the initial source is irro- 
tational. Taking Fourier transform of both sides in (3.26), we obtain that 

u B (k;x )= a,(r,xo,k)^ l (k)(Ao(xo)rb rl (xo,k)-b,(xo,k))rb 3 (x ,k) (3.28) 

j,n£j 



where we have defined 



S(k) = ^£„(k)b„(x ,k). 



(3.29) 



n£j 



Irrotationality of the initial source means that Si and S2 identically vanish, or equiv- 
alent^ that 



S(x) 



/V</>(x) 



(3.30) 

-bzp and by 



for some pressure p(x) and potential </>(x). Remarking that rb± 
irrotationality that (-Ao(xo)S(k) • bi,2(k)) = 0, we use (3.23) to recast (3.28) as 

u B (k; x ) = a_ (T, x , k)S+ (k)b+ (x , k) + a+ (T, x , k) S- (k)b_ (x , k). (3.31) 
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Decomposing the source S(x) as 

S(x) = S+(x) + S_(x), such that S±(k) = 1 S , ± (k)b ± (x , k), 
the back-propagated signal takes the form 

u B (|;x ) = (a_(T,x ,-) * S+(-))(£) + («+(T,x , •) * S_(-))(€) (3.32) 



where a± is the Fourier of a± in k. This form is much more tractable than (3.26) 



( 3.27 ). It is also almost as general. Indeed, rotational modes do not propagate 
in the high frequency regime. Therefore, they are exactly back-propagated when 
x(xo) = 1 and /(x) = <5(x), and not back-propagated at all when x( x o) = 0. All 
the refocusing properties are thus captured by the amplitudes a±(T, Xo,k). Their 
evolution equation characterizes how waves propagate in the medium and their initial 
conditions characterize the recording array. 

3.6. Homogeneous Media. In homogeneous media with c(x) = cq the ampli- 
tudes a± (T, x, k) satisfy the free transport equation Jl4|, Q 

-J± ± c k ■ V x a± = (3.33) 
with initial data a±(0, x, k) = |x(x)| 2 /(k) as in ( |3.25 ). They are therefore given by 



a±(i,x ,k) = | x (x oT c kt)| 2 /(k). (3.34) 

These amplitudes become more and more singular in k as time grows since their 
gradient in k grows linearly with time. The corresponding kernel F = Fh decays 
therefore more slowly in £ as time grows. This implies that the quality of the refocus- 
ing degrades with time. For sufficiently large times, all the energy has left the domain 
fl (assumed to be bounded), and the coefficients a±(t, xo,k) vanish. Therefore the 
back-propagated signal u B (£; xq) also vanishes, which means that the re is no refocus- 



ing at all. The same conclusions could also be drawn by analyzing (3A) directly in 
a homo gen eous medium. This is the situation in the numerical experiment presented 
in Fig. p.2[ in a homogeneous medium, the back-propagated signal would vanish. 

3.7. Heterogeneous Media and Radiative Transport Regime. The results 
of the preceding sections show how the back-propagated signal u B (£;xo) is related 
to the propagating modes a±(T, xo,k) of the Wigner matrix W(T, xo,k). The form 
assumed by the modes a±(T, Xo,k), and in particular their smoothness in k, will 
depend on the hypotheses we make on the underlying medium; i.e., on the density 
p(x) and compressibility k(x) that appear in the matrix A(x). We have seen that 
partial measurements in homogeneous media yield poor refocusing properties. We 
now show that refocusing is much better in random media. 

We consider here the radiative transport regime, also known as weak coupling 
limit. There, the fluctuations in the physical parameters are weak and vary on a scale 
comparable to the scale of the initial source. Density and compressibility assume the 
form 

x x 

/?(x) = p + v / ipi(-) and k(x) = k + VeKi(-). (3.35) 

e e 

The functions p\ and K\ arc assumed to be mean-zero spatially homogeneous pro- 
cesses. The average (with respect to realizations of the medium) of the propagating 
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amplitudes a±, denoted by a±, satisfy in the high frequency limit E->0a radiative 
transfer equation (RTE), which is a linear Boltzmann equation of the form 

da± 



dt 



± Cok • V x a± 



<x(k,p)(a ± (i,x,p) - a±(t, x, k))5(co(|k| - |p|))dp 



(3.36) 



a±(0,x,k) = | X (x)| 2 /(k). 



The scattering coefficient cr(k, p) depends on the power spectra of p\ and K\. We refer 
to |24| for the details of the derivation and explicit form of cr(k, p). The above result 
remains formal for the wave equation and requires to average over the realizations 
of the random medium although this is not necessary in the physical and numerical 
time reversal experiments. A rigorous proof of the derivation of the linear Boltzmann 
equation (which also requires to average over realizations) has only been obtained 
for the Schrodinger equation; see || g6). Neve rthel ess, t he ab ove result formally 
characterizes the filter F(T, £; xp) introduced in ( 3.27 ) and ( |3.32[ ). 

The transport equation (3.36) has a smoothing effect best seen in its integral 
formulation. Let us define the total sca tterin g coefficient S(k) = J R3 cr(k, p)<5(co(|k| — 
|p| ))dp. Then the transport equation ( 3.36| ) may be rewritten as 



a±(t, x, k) = a±(0, x =p c kt, k)e~ s(k)t 

<r(k, |k|p)a±(s,x =f c (t - s)k, |k|p)e~ 



(3.37) 



+ 



Co 



ds 



E(k)(i-«) 



dfi(p). 



Here p = p/|p| is the unit vector in di recti on of p and d£l(p) is the surface element 
on the sphere S 2 . The first term in ( |3.37| ) is the ballistic part that undergoes no 
scattering. It has no sm ooth ing effect, and, moreover, if a(0, x, k) is not smooth in x, 
as may be the case for (|3.25 ), the discontinuities in x translate into discontinuities in 
k at latter times as in ( p.34| ) in a homogeneous medium. However, in contrast to the 
homogeneous medium case, the ballistic term decays exponentially in time, and does 
not affect the refocused signal for sufficiently long times t 3> 1/E. The second term 
in (3.37) exhibits a smoothing effect. Namely the operator Cg defined by 



Cg{t, x, k) = M f ds f cr(k, \k\p)g(s, x t c (t - s)k, |k|p)e- s ( k )(*- s )^(p) 
c o Jo Js 2 

is regularizing, in the sense that the function g = Cg has at least 1/2-more derivatives 
than g (in some Sobolev scale). The precise formulation of this smoothing property 
is given by the averaging lemmas ]l5| , |22| and will not be dwelt upon here. Iterating 
(3.37) n times we obtain 

■ • + <4 (t, x, k) + £ n+1 a± (t, x, k) . (3.38) 



a± (t, x, k) = a° (t, x, k) + o± (t, x, k) + ■ 
The terms a^. , . . . , a± are given by 

4(t,x,k) = a±(0,x=Fc o kt,k)e" s(k) ', 



a±(t, x, k) 



(*,x,k) 



They describe, respectively, the contributions from waves that do not scatter, scatter 
once, twice, .... It is straightforward to verify that all these terms decay exponen- 
tially in time and are negligible for times t ^ 1/E. The last term in ( 3.38| ) has 
at least n/2 more derivatives than the initial data ao, or the solution (3.34) of the 
homogeneous transport equation. This leads to a faster decay in £ of the Fourier 
transforms a±(T, Xo,£) of a±(T, Xo,k) in k. This gives a qualitative explanation as 
to why rcfocusing is better in heterogeneous media than in homo geneo us media. A 
more quantitative answer requires to solve the transport equation ( 3.36 ). 
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3.8. Diffusion Regime. It is known for times t much longer than the scattering 
mean free time r sc = 1/S and distances of propagation L very large compared to l sc = 
cqt sc that solutions to the radiative transport equation ( 3.36| ) can be approximated by 
solutions to a diffusion equation, provided that c(x) = cq is independent of x [2^] . 
More precisely we let S — l sc /L C 1 be a small parameter and rescale time and 
space variables as t — -> t/S 2 and x — > x/<5. In this limit, wave direction is completely 
randomized so that 



a + (t, x, k) ?s a_ (t, x, k) ps a(t, x, |k| 



where a solves 



da(t, x, |k| 



£>(|k|)A x o(t,x,|k|) = 0, 



a(0,x,|k|) = | X (x)p 



47r|k|' 



/(q)<J(|q| - |k|)dq. 



(3.39) 



The diffusion coefficient £)(|k|) may be expressed explicitly in terms of the scattering 
coefficient <r(k, p) and hence related to the power spectra of pi and K\. We refer to 
[PH for the details. For instance, let us assume for simplicity that the density is not 
fluctuating, pi = 0, and that the compressibility fluctuations are delta-correlated, so 
that E{Ki(p)«i(q)} = KqRq5(p + q). Then we have 



cr(k,p) 



nc 2 \k\ 2 R 



and 



Eflkl) = 27r 2 c |k| 4 i? 



CO 



3E(|k|) 67r 2 |k| 4 i? 



(3.40) 



(3.41) 



Let us assume that there are no initial rotational modes, so that the source S(x) 
is decomposed as in ( |3~30|) . Using ( |3~3i"l ), we obtain that 



u B (k;x ) =o(r,xo,|k|)S(k). 



(3.42) 



When /(x) is isotropic so that /(k) = /(|k|), and the diffusion coefficient is given by 
( |3.41 ), the solution of (3.3E) takes the form 



o(T,xo,|k|)=/(|k|) 



37r|k| 4 i? \ 3 / 2 
2c T 



exp 



3^ 2 |k| 4 i? |x -y|' 
2c T 



\x(y)\ 2 dy- 



(3.43) 



When /(x) = S(x), and = K 3 , so that x(x) = 1, we retrieve a(T, xo,k) = 1, hence 
the refocusing is perfect. When only partial measurement is available, the above 
formula indicates how the frequencies of the initial pulse arc filtered by the one-step 
time reversal process. Notice that both the low and high frequencies are damped. 
The reason is that low frequencies scatter little with the underlying medium so that 
it takes a long time for them to be randomized. High frequencies strongly scatter 
with the underlying medium and consequently propagate little so that the signal that 
reaches the recording array fl is small unless recorders are also located at the source 
point: xq € f2. In the latter case they are very well measured and back-p ropa gated 
although this situation is not the most interesting physically. Expression ( 3.4S ) may 
be generalized to other power spectra of medium fluctuations in a straightforward 
manner using the formula for the diffusion coefficient in ||24| . 
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3.9. Numerical Results. The numerical results in Fig. 2.2 show that some 
signal refocuses at the location of the initial source after the time reversal procedure. 
Based on the above theory however, we do not expect the refocused signal to have 
exactly the same shape as the original one. Since the location of the initial source 
belongs to the recording array (x(xo) = 1) in our simulations, we expect from our 
theory that high frequencies will rcfocus well but that low frequencies will not. This is 



Fig. 



2iq 

2. 



3.1. Zoom of the initial source and the refocused signal for the numerical experiment of 



confirmed by the numerical results in Fig. 3.1 , where a zoom in the vicinity of x = 
of the initial source and refocused signal are represented. Notice that the numerical 
simulations are presented here only to help in the understanding of the refocusing 
theory and do not aim at reproducing the theory in a quantitative manner. The 
random fluctuations are quite strong in our numerical simulations and it is unlikely 
that the diffusive regime may be valid. The refocused signal on the right figure looks 
however like a high-pass filter of the signal on the left figure, as expected from theory. 

4. Refocusing of Classical Waves. The theory presented in section|| provides 
a quantitative explanation for the results observed in time reversal physical and nu- 
merical experiments. However, the time reversal procedure is by no means necessary 
to obtain refocusing. Time reversal is associated with the specific choice ( |3.2|) for the 
matrix T in the preceding section, which reverses the direction of the acoustic veloc- 
ity and keeps pressure unchanged. Other choices for T are however possible. When 
nothing is done at time T, i.e., when we choose r = /, no refocusing occurs as one 
might expect. It turns out that T = I is more or less the only choice of a matrix that 
prevents some sort of refocusing. Section 4.1 presents the theory of refocusing for 
acoustic waves, which is corroborated by numerical results presented in section 4.2. 
Sections 4.3 and 4.4 generalize the theory to other linear hyperbolic systems. 



4.1. General Refocusing of Acoustic Waves. In one-step time reversal, the 
action of th e "in telligent" array is captured by the choice of the signal processing 
matrix T in (3.3). Time reversal is characterized by T given in (3^2). A passive array 
is characterized by T = I. This section analyzes the role of other choices for T, which 
we let depend on the receiver location so that each receiver may perform its own kind 
of signal processing. 
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The signal after time reversal is still given by ( |3.3| ), where r(y') is now arbi- 
trary. At time 2T, after back-propagation, we are free to multiply the signal by an 
arbitrary invertible matrix to analyze the signal. It is convenient to multiply the back- 
propagated signal by the matrix Tq = Diag(— 1, —1, —1, 1) as in classical time reversal. 
The reconstruction formula (3.5) in the localized source limit is then replaced by 

u s (£;x )= / roG(r,xo+^;y)r(y')G(T,y';x +ez)x(y,y')S( Z )dydy'd Z (4.1) 

with x(y, y') defined by (|3.6|). To generalize the results of section ||, we need to define 
an appropriate adjoint Green's matrix G*. As before, this will allow us to remove 
the matrix T between the two Green's matrices in (fOl) and to interchange the order 
of points in the second Green's matrix. We define the new adjoint Green's function 
G* (t, x; y) as the solution to 



9G*(i,x;y) 



A(x) 



dG*(t,x;y) 



= 



dt dxi 
G* (0, x; y) = r(x)r A' 1 (x)<5(x - y) 



Following the steps of section 3.3, we show that 

G,(t,x,y) =r(y)G(i,y;x)A- 1 (x)r c 



(4.2) 



(4.3) 



The only modification compared to the corresponding derivation of (3.8) is to multiply 

(3.9) on the left by T(x) and on the right by Tq so that T(y) appears on the left in 

(3.10) . The re-transmitted signal may now be recast as 

u B (£;x )= f rfydy'dzr G(T,x + ^;y)G4r,xo+£z;y')r- 1 (4.4) 

x ^(x + ez)x(y,y')S(z). 

Therefore the only modification in the expression for the re-transmitted signal com- 
pared to the time reversed signal (3.12) is in the initial data for (4.2), which is the 
only place where the matrix T(x) appears. 

The analysis in sections |3 . 3| - [3 . 7| requires only minor changes, which we now outline. 
The back-propagated signal may still be expressed in term of the Wigner distribution 
(compare to (3.17)) 



u B (£; xo) = / e^-^T a W £ (T, x + ^±1, k)r A(x + ez)S(z) 



dzdk 

(2^)3' 



(4.5) 



The Wigner distribution is defined as before by (3.15) and ( 3.1 6| ) . The function Q is 
defined as before as the solution of ( f^3| ) with initial data Q(0,a;;q) = x(x)e lq x / £ /, 
while Q* solves with the initial data Q»(0,x;q) = r(x)r A _1 (x)x(x)e- iq ' x / E . 
The initial Wigner distribution is now given by 

W(0,x,k) = | X (x)| 2 r(x)r ^- 1 (x)/(k). (4.6) 

Lemma 3.1 and Proposition 3.2 also hold, and we obtain the analog of ( 3.19| ) 



u(£;x ) 



^<i-')T W(T, x , k)r A (x )S(z)dzdk. 



(4.7) 
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The limit Wigner distribution W(T, x ,k) admits the mo de d ecomposition ( 3.24 ) as 
before. If we assume that the source S(x) has the form ( 3.3C ) so that no rotational 
modes are present initially, we recover the refocalization formula ( [3.31 ) : 



u s (k;x ) = a_(T, x ,k)5 + (k)b+ (x ,k) +a + (T, x ,k)S , _(k)b_(x ,k). (4.8) 



The initial conditions for the amplitudes a± are replaced by 

o± (0, x, k) = Tr [A (x) W(Q, x, k) A Q (x)b± (x , k)b^ (x , k) 
= | X (x)| 2 /»(^o(x)r(x)b T (x,k) • b±(x,k)). 



(4.9) 



Observe that when T(x) = To, we get back the results of section 3.7. When the 
signal is not changed at the array, so that T = /, the coefficients a±(0,x, k) = by 
orthogonality (3.23) of the eigenvectors hj. We thus obtain that no refocusing occurs 



when the "intelligent" array is replaced by a passive array, as expected physically. 

Another interesting example is when only pressure p is measured, so that the 
matrix T = Diag(0, 0, 0, 1). Then the initial data is 

a± (0,x,k) = i| X (x)| 2 />), 



which differs by a factor 1/2 from the full time reversal case (3.25). Therefore the 
re- transmitted signal u B also differs only by a factor 1/2 from the latter case, and 
the quality of refocusing as well as the shape of the re-propagated signal are exactly 
the same. The same observation applies to the measurement and reversal of the 
acoustic velocity only, which corresponds to the matrix T = Diag(— 1, —1,-1,0). The 
factor 1/2 comes from the fact that only the potential energy or the kinetic energy 
is measured in the first and second cases, respectively. For high frequency acoustic 
waves, the potential and kinetic energies are equal, hence the factor 1/2. We can also 
verify that when only the first component of the velocity field is measured so that 
r = Diag(— 1, 0, 0, 0), the initial data is 



a±(0,x,k) = |x(x)| 2 /(k) 



2|kl : 



(4.10) 



As in the time reversal setting of section [| the quality of the refocusing is related 
to the smoothness of the amplitudes a± in k. In a homogeneous medium they satisfy 
the free transport equation ( |3. 33 ) , and are given by 



a±(i,x,k)-| X (x- Co kt)| 2 /(k) 

x (Aq(x — c ki)r(x — coki)b T (x — cok£, k) • b±(x — coki, k)). 

Once again, we observe that in a uniform medium a± become less regular in k as time 
grows, thus refocusing is poor. 



The considerations of section 3/7 show that in the radiative tran spor t regime the 
amplitudes a± become smoother in k also with ini tial data given by (4.9). This leads 
to a better refoc using as explained in section 3J: . Let us assume that the diffusion 
regime of section 3J3 is valid and that the kernel / is isotropic /(k) = /(|k|). This 
requires in particular that An(x) be independent of x. We obtain that a±(T, Xo, k) = 
a(T, xo, |k|), thus the refocusing formula ( |4.8| ) reduces to 



u B (k;x ) =a(T, x , |k|)S(k). 



(4.11) 
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The difference with the case treated in section 3.8 is that a(T, x, |k|) solves the diffu- 
sion equation (3.39) with new initial conditions given by 



a(0,x,|k|) = M^ / 



47r|k| 2 

lx(x)| 2 
47r|k| 2 



/(|q|)( J 4 r(x)b + (q) • b_(q))S(|q| - |k|)dq. 



When only the first component of the velocity field is measured, as in fl41C| ), the 
initial data for a is 



a(0,x,|k|) = -| X (x)| 2 /(|k|). 

Therefore even time reversing only one component of the acoustic velocity field pro- 
duces a re-propagated signal that is equal to the full re-propagated field up to a 
constant factor. 



More generally, we deduce from (4.12) that a detector at x will contribute some 
refocusing for waves with wavenumber |k| provided that 



/(|k|q)(A r(x)b T (q) • b±(q))dfi(q) ^ 0. 



When /(x) = /(|x|) is radial, this property becomes independent of the wavenumber 
|k| and reduces to j s2 (A r(x)b T (q) • b ± (q))cif2(q) ^ 0. 

4.2. Numerical Results. Let us come back to the numerical results presented 
in Fig. [T^ and 3.1. We now consider two different processings at the recording 
array. The first array is passive, corresponding to T = I, and the second array only 
measures pressure so that T = Diag(0, 0, 0, 1). The zoom in the vicinity of xo = of 
the "refocused" signals is given in Fig. 4.1. The left figure shows no refocusing, in 



Fig. 4.1. Zoom of the refocused signals for the numerical experiment of Fig. 2.i with processing 
T = I (left), with a maximal amplitude of roughly 4 10 — 3 and F = Diag(0, 0, 0, 1) (right), with a 
maximal amplitude of roughly 0.035. 



accordance with physical intuition and theory. The right figure shows that refocusing 
indeed occurs when only pressure in recorded (and its time derivative is set to in 
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the solution of the wave equation presented in the appendix), 
refocuscd signal is roughly one half the one obtained in Fig. 
theory. 



Notice also that the 
pTl] as predicted by 



4.3. Refocusing of Other Classical Waves. The preceding sections deal with 
the refocusing of acoustic waves. The theory can however be extended to more com- 
plicated linear hyperbolic systems of the form with A(x) a positive definite 
matrix, D- 7 symmetric matrices, and u <E C m . These i nclu de electromagnetic and 
elastic waves. Their explicit representation in the for m (|2.3|) and expressions for the 
matrices A(x) and Z? J in these cases may be found in 
equations 




For instance, the Maxwell 



<9E 
~dt 



1 



e(x) 



curl H 



rl E 



may be written in the form ( |2.3| ) with u = (E, H) £ C 6 and the matrix A(x) = 
Diag(e(x), e(x), e(x), /n(x), /x(x), /x(x)). Here e(x) is the dielectric constant (not to be 
confused with the small parameter e), and /i(x) is the magnetic permeability. The 
6x6 dispersion matrix L(x, k) for the Maxwell equations is given by 



L(x,k) = - 








-fe 3 //i(x) 







fc 3 //i(x) 


-fci//i(x) 







-fc 2 //i(x) 
fc x //i(x) 






fc 3 /e(x) 
-fc 2 /e(x) 






-*3/e(x) 


fci/e(x) 






fc 2 /e(x) \ 
-fex/e(x) 









Generalization of our results for acoustic waves to such general systems is quite 
straightforward so we concentrate only on the modifications that need be made. The 
time reversal procedure is exactly the same as before: a signal propagates from a 
localized source, is recorded, processed as in (3.3) with a general matrix r(y'), and 
re-emitted into the medium. The re-transmitted signal is given by (4.1). Further- 
more, the equation for the adjoint Green's matrix ( |4.2| ), the definition of the Wigncr 
transform in section 3.4, and the expression (4.7) for the re-propagated signal still 
hold. 

The analysis of the re-propagated signal is reduced to the study of the Wigncr 
distribution, which is now modified. The mode decomposition need be generalized. 
We recall that 

L(x,k)=A - 1 (x)fc,^ 

is the m x m dispersion matrix associated with the hyperbolic system ( |2.3| ). Since 
L(x, k) is symmetric with respect to the inner product (u, v)^,, = (Aq\i ■ v), its 
eigenvalues are real and its eigenvectors form a basis. We assume the existence of a 



time reversal matrix To such that (3.11) holds with T = To and such that Tq = /. For 
example, for electromagnetic waves Tq = Diag(l, 1,1,-1,-1,-1). Then the spectrum 
of L is symmetric about zero and the eigenvalues ±cj q have the same multiplicity. We 
assume in addition that L is isotropic so that its eigenvalues have the form oj± (x, k) = 
±c"(x)|k|, where c Q (x) is the speed of mode a. We denote by r a their respective 
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multiplicities, assumed to be independent of x and k for k ^ 0. The matrix L has a 
basis of eigenvectors b ^ (x, k) such that 

L(x,k)b^'(x,k) = ±^(x,k)b^'(x,k), j = l,...,r a , 

and h± J form an orthonormal set with respect to the inner product (,}a - The 
different uj a correspond to different types of waves (modes) . Various indices 1 < j < r a 
refer to different polarizations of a given mode. The eigenvectors h°t_' J and b"' J are 
related by 

r b^(x,k) = b^'(x,k), r b^'(x,k) = b^'(x,k). (4.13) 



Proposition 3.2 is then generalized as follows |14| |24|: 

Proposition 4.1. There exist scalar functions a^ : ' :,m (i, x, k) such that 

W(t,x,k)= a± Jm (^x,k)b^'(x,k)®bJ m (x,k). (4.14) 

Here the sum runs over all possible values of ±, a, and 1 < j,m < r a . 

The main content of this proposition is again that the cross terms b ^ J (x, k) <E> 

b^' m (x, k) do not contribute, as well as the terms \x± 3 (x, k)(g)b" '™(x, k) when a ^= a' . 
This is because modes propagating with different speeds do not interfere constructively 
in the high frequency limit. 

We m ay n ow insert expression ( 4.14 ) into (4/7) and obtain the following general- 
ization of (P~8|) 

u B (k;x )= £ [ a ^'(T,x ,k)^ J (x ,k)b^ m (x ,k) (4.15) 

ct,j,m 

+ a«> mj (T, xo,k) 1 9^(x 0) k)b^ m (x ,k)" , 

where S , ^ J (k) = (A(xo)S(k) • b± J (xo, k)). This formula tells us that only the modes 
that are present in the initial source (S^ 1 (k) ^= 0) will be present in the back- 
propagated signal but possibly with a different polarization, that is, j ^ m. 
The initial conditions for the modes a± jm are given by 

a^ m (0,x,k) = | X (x)| 2 /(k)(A(x)r(x)b^ m (x,k) • b^'(x,k)), (4.16) 

which generalizes (p~9|). When T(x) = /, we again obtain that a± jm (0, x, k) = 0, i.e., 
there is no refocusing as physically expected. When T(x) = To, we have for all a that 

a^ m (0,x,k) = | X (x)| 2 /»<5, m . 

In a uniform me dium the amplitudes a^. jm satisfy an uncoupled system of free trans- 
port equations ( 3.33 ): 

Pt„oc,jm 

± c Q k • V x a^ m = 0, (4.17) 



dt 

which have no smoothing effect, and hence refocusing in a homogeneous medium is 
still poor. When /(x) = <5(x) and Q, = R 3 , so that x(x) = 1, we still have that 
a °^> m (T', xo, k) = S j m a nd refocusing is again perfect, that is, u B (£;xo) = S(£), as 
may be seen from (4.15). 
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4.4. The diffusive regime. The radiative transport regime holds when the 
matrices A(x) have the form 



as in (3.35). Then the r a x r a coherence matrices w± with entries w± 



Jin 



satisfy a system of matrix- valued radiative transport equations (see E4j for the details) 



similar to (3.36). The matrix transport equations simplify considerably in the diffusive 



regime, such as the one considered in section when waves propagate over large 
distances and long times. We assume for simplicity that Aq = Aq(x) and T = T(x) 



are independent of x. Polarization is lost in this regime, that is, 



l (i,x,k) = for 



j ^ m and wave energy is equidistributed over all directions. This implies that 



(t,x,k) 



a, 77 



(i,x, k) = a«(i,x, |k| 



so that a a is independent of j = 1, . . . , r a and of the direction k = k/|k|. Further- 
more, because of multiple scattering, a universal equipartition regime takes place so 
that 



O a (t,Xo, |k|) = 0(t,Xo,CQ,|k|), 



(4.18) 



where <j)(t, x, ui) solves a diffusion equation in x like ( |3~39l) (sec ||). The diffusion 
coefficient D(uS) may be expressed explicitly in terms of the power spectra of the 
medium fluctuations |24|j . Using ( 4.16| ) and ( |l.l§| ), we obtain when / is isotropic the 
following initial data for the function <p 



0(0, x,w) 



ilx(x)l 2 / ^ E /(^)(^orb-(k),b-(k)) d o ( k )) 



(4.19) 



where \a\ is the number of non-vanishing eigenvalues of L(x, k), and dfl(k) is the 
Lebesgue measure on the unit sphere S* 2 . 

Let us assume that non-propagating modes are absent in the initial source S(x), 
with the subscript zero referring to modes corresponding to ujq = 0. 



that is, S^(k 



Then ( 4.15 ) becomes 



u(k;x ) =£>CZ> ,c a |k|) [S^(k)b^'(x ,k) + S^(k)b^'(x ,k) 



(4.20) 



This is an explicit expression for the re-propagated signal in the diffusive regime, 
where <j> solves the diffusion equation ( 3.39| ) with initial conditions ( 4.19 ). 



5. Conclusions. This paper presents a theory that quantitatively describes the 
refocusing phenomena in time reversal acoustics as well as for more general processings 
of other classical waves. We show that the back-propagated signal may be expressed 



as the convolution (1.1) of the original source S with a filter F. The quality of the 
refocusing is therefore determined by the spatial decay of the kernel F. For acoustic 



waves, the explicit expression (3.27) relates F to the Wigner distribution of certain 
solutions of the wave equation. The decay of F is related to the smoothness in the 
phase space of the amplitudes aj(t, x, k) defined in Proposition 3.3. The latter satisfy 
a free transport equation in homogeneous media, which sharpens the gradients of cij 
and leads to poor refocusing. In contrast, the amplitudes dj satisfy the radiative 
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transport equation ( 3.36| ) in heterogeneous media, which has a smoothing effect. This 



leads to a rapid spatial decay of the filter F and a better refocusing. For l onge r t imes , 
satisfies a diffusion equation. This allows for an explicit expression ( 3.42] )-( 3.43| ) 



of the time reversed signal. The same theory holds for more general waves and more 
general processing procedures at the recording array, which allows us to describe the 
refocusing of electromagnetic waves when only one component of the electric field is 
measured, for instance. 

Appendix. This appendix presents the details of the numerical simulation of 
( |2.9| ). We assume that p is constant and that only k(x) fluctuates. We can therefore 
recast ( pljj ) as 

- c 2 (x)A P = 0. 

The above wave equation is discretized using a second-order scheme (three point 
stencil in every variable) both in time and space. The resolution in time is explicit 
and time reversible, i.e., the equation that yields p(t n+ i) fromp(i„_i) andp(t„) can be 
used to retrieve p(t n -i) exactly from p(t n ) and p(t n+ i). We write c 2 (x) = Cq + c\ (x). 
The average velocity is Cq = 1. The random part c 2 has been constructed as follows. 
Let 2N x 2N be the number of spatial grid points and c 2 .„ m be the value of cf at the 
grid point (n, m). The values c 2 . 2n 2m have been chosen independently and uniformly 
on (— r, r) with r < 1/2. The value of c 2 is then set constant on four adjacent pixels by 
enforcing that cf. 2n _ 1)2m = c 2 . 2 „_ lj2m _ 1 = c?. an>2m _i = c 2 ;2n ,2m for 1 < n,m < N. 
In all simulations, we have N = 200, which generates a grid of 400 2 = 1.6 10 4 points. 
The time step has been chosen so that the CFL condition St < min c(x) / (2N) is 

ensured. 
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